!======================================================================
 module gazlib
!---------------------------
           contains
!  SUBROUTINE MULTY(EXT,RES,NN,M,N,EPS0)
!  SUBROUTINE DADAPT(F,A,B,NSEG,RELTOL,ABSTOL,RES,ERR)
!  real(8) FUNCTION DGAUSS(F,A,B,EPS)
!  Subroutine GAUSS(F,A,B,EPS,H)
!  SUBROUTINE QUANC8(FUN,A,B,ABSERR,RELERR,RESULT,ERREST,NOFUN,FLAG)
!  real(8) Function FLin(x,xa,xb,ya,yb)
!  real(8) FUNCTION DIVDIF(F,X,N,Z,M)
!  SUBROUTINE EXTVEC(X,Y,F,N,M,L,XX,YY,FF)
!  Subroutine FindPos(N, XX, x, i)
!  SUBROUTINE DDEQMR(N,XA,XZ,Y,H0,EPS,SUB,W)
!  real(8) Function ERF(X)
!  real(8) FUNCTION BESSK0(X)
!  real(8) FUNCTION BESSK1(X)
!  real(8) FUNCTION BESSK2(X)
!  real(8) FUNCTION BESSI0(X)
!  real(8) FUNCTION BESSI0(X)

!---------------------------------

!==========================================
! 05.12.12, LNGS
Subroutine cprb(x,y,c)
   implicit none
   real(8),intent( in),dimension(3):: x, y ! input values of [x] and [y]
                                           ! derivatives of parab at xv
   real(8),intent(out),dimension(3):: c    ! array of coefficients

   real(8):: x12, x23, x13
   
   x12  =   x(1) - x(2)
   x23  =   x(2) - x(3)
   x13  =   x(1) - x(3)

   c(1) =   x(2)*y(1) - x(3)*y(1) - x(1)*y(2)      &
          + x(3)*y(2) + x(1)*y(3) - x(2)*y(3)
   
   c(2) = - x(2)**2*y(1) + x(3)**2*y(1)            &
          + x(1)**2*y(2) - x(3)**2*y(2)            &
          - x(1)**2*y(3) + x(2)**2*y(3)            
	     
   c(3) =   x(2)**2*x(3)*y(1) - x(2)*x(3)**2*y(1)  &
          - x(1)**2*x(3)*y(2) + x(1)*x(3)**2*y(2)  &
		  + x(1)**2*x(2)*y(3) - x(1)*x(2)**2*y(3) 
		  
   c    = c/(x12*x13*x23)

end Subroutine cprb
!==========================================

!==========================================
! 05.12.12, LNGS
Subroutine 	parab(n,xx,yy,x,y,ys,yss)
   implicit none
   integer,intent(in):: n
   real(8),intent(in),dimension(n):: xx, yy
   real(8),intent(in):: x
   real(8),intent(out):: y, ys, yss     ! output value, first and second 
   real(8),dimension(3):: x3, y3
   real(8),dimension(3):: c             ! array of coefficients
   integer:: i, ibeg

   call FindPos(n,xx,x,ibeg)
   if(ibeg.eq.n-1) ibeg = ibeg-1
   if(ibeg.eq.n) ibeg = ibeg-2

   do i = 1,3
      x3(i) = xx(ibeg+i-1)
	  y3(i) = yy(ibeg+i-1)
   enddo ! i
   
   call cprb(x3,y3,c)
!   print*,' c=',c
   y   = (c(1)*x + c(2))*x + c(3)
   yss = 2d0*c(1)
   ys  = yss*x + c(2)

end Subroutine parab
!==========================================

!***************************************************************
!                   INTEGRATION
!***************************************************************

!==========================================
! 28.11.12, LNGS
! rewritten for FORTRAN 90

SUBROUTINE MULTY(EXT,RES,NN,M,N,EPS0)
   implicit none
   integer,intent(in):: M, N
   integer,intent(in),dimension(N):: NN
   real(8),dimension(M):: RES
   real(8),intent(in):: EPS0
   real(8):: EPS, UK, AA, BB, DW
 
   INTEGER P, Q
   integer,dimension(4):: H, K, S  
   integer:: L, I, I1, J, JJ, K1, K2, E, E1

   real(8),DIMENSION( 4):: A, D, X, U, W
   real(8),DIMENSION(10):: Y, V
   real(8),DIMENSION(4,10):: F, R
   
   DATA U/0.8611363d0, 0.3399810d0, 0d0, 0d0/   &
      , W/0.3478548d0, 0.6521452d0, 0d0, 0d0/   &
      , P/4/
!- - - - - - - - - - - -
   EPS=EPS0
   DO JJ=1,N  
      S(JJ)=MAX0(1,NN(JJ))
   enddo  ! JJ

   DO I=1,M
      V(I)=0d0
   enddo
   
   L=0
   DO JJ=1,N
14    IF(L.EQ.1) S(JJ)=2*S(JJ)
      L = 1
      J = 0
      GOTO 3
!- - - - - - - - - - - -
1     K(J) = 1
2     K1 = K(J)
      K2 = P+1-K1
      UK = -U(K1)
      IF(K1.GT.K2) UK = U(K2)
      X(J) = A(J) + (2d0*H(J) - 1d0 + UK)*D(J)
      IF(J.EQ.N) GOTO 5
3     Q = J+1
      DO J=Q,N
         CALL EXT(J,X,AA,BB,Y)
         A(J) = AA
         D(J) = (BB-AA)/(2*S(N+1-J))
         X(J) = AA + (1d0 - U(1))*D(J)
         H(J) = 1
         K(J) = 1
         DO I=1,M
            R(J,I) = 0d0
            F(J,I) = Y(I)
         enddo ! I
      enddo ! J  
      J = N
5     CALL EXT(N+1,X,AA,BB,Y)
      K1 = MIN0(K(J),P +1 - K(J))
      DO I=1,M
         R(N,I) = R(N,I) + Y(I)*W(K1)
      enddo ! I
7     K(J) = K(J) + 1
      IF (K(J).LE.P) GOTO 2
      H(J) = H(J) + 1
      IF(H(J).LE.S(N+1-J)) GOTO 1
      J=J-1
      IF(J.EQ.0) GOTO 9
      K1=MIN0(K(J),P+1-K(J))
      DW = D(J+1)*W(K1)
      DO I=1,M
         R(J,I) = R(J,I) + R(J+1,I)*F(J+1,I)*DW
      enddo ! I
      GOTO 7
9     DO I=1,M
         Y  (I) = S(I)
         RES(I) = R(1,I)*F(1,I)*D(1)
      enddo ! I
!     - - - - - - - - - - - -
      CALL EXT(N+2,Y,EPS,EPS,RES)
      IF(EPS.LE.0d0) RETURN
      DO I=1,M
         E  = DABS(RES(I)-V(I))
         E1 = DABS(EPS0*RES(I))
         IF(E.LE.E1) EXIT
         DO I1=1,M
            V(I1) = RES(I1)
         enddo ! I1
         GOTO 14
      enddo ! I
      S(JJ) = S(JJ)/2
   enddo ! JJ
   RETURN
END SUBROUTINE MULTY

!**************************************************************!
      SUBROUTINE DADAPT(F,A,B,NSEG,RELTOL,ABSTOL,RES,ERR)
      IMPLICIT none
      real(8),intent(in):: A, B, RELTOL, ABSTOL 
      real(8),intent(out):: RES, ERR
      integer:: I, NSEG, NSEGD, ITER, IBIG, NTER=0
      integer, parameter:: NDIM=100
      real(8), parameter:: R1 = 1d0, HF = R1/2d0

      real(8), dimension(NDIM):: XLO,XHI,TVAL,TERS
      real(8):: F, TVALS, TERSS, TE, ROOT, XLOB, XNEW &
              , XHIB, BIN, BIGE
      external F

      SAVE XLO,XHI,TVAL,TERS,NTER

!     RES = Estimated Integral of F from A to B,
!     ERR = Estimated absolute error on RES.
!     NSEG  specifies how the adaptation is to be done:
!        =0   means use previous binning,
!        =1   means fully automatic, adapt until tolerance attained.
!        =n>1 means first split interval into n equal segments,
!             then adapt as necessary to attain tolerance.
!     The specified tolerances are:
!            relative: RELTOL ;  absolute: ABSTOL.
!        It stops when one OR the other is satisfied, or number of
!        segments exceeds NDIM.  Either TOLA or TOLR (but not both!)
!        can be set to zero, in which case only the other is used.
  
  
      IF(NSEG .LE. 0)  THEN
         IF(NTER .EQ. 0) THEN
            NSEGD=1
            GO TO 2
         ENDIF
         TVALS=0d0
         TERSS=0d0
         DO I = 1,NTER
            CALL DGS56P(F,XLO(I),XHI(I),TVAL(I),TE)
            TERS(I)=TE**2
            TVALS=TVALS+TVAL(I)
            TERSS=TERSS+TERS(I)
         ENDDO
         ROOT= DSQRT(2d0*TERSS)
         GO TO 9
      ENDIF

      NSEGD=MIN(NSEG,NDIM)
    2 XHIB=A
      BIN=(B-A)/NSEGD
      DO I = 1,NSEGD
         XLO(I)=XHIB
         XLOB=XLO(I)
         XHI(I)=XHIB+BIN
         IF(I .EQ. NSEGD) XHI(I)=B
         XHIB=XHI(I)
         CALL DGS56P(F,XLOB,XHIB,TVAL(I),TE)
         TERS(I)=TE**2
      ENDDO

      NTER=NSEGD
      DO ITER = 1,NDIM
         TVALS=TVAL(1)
         TERSS=TERS(1)
         DO I = 2,NTER
            TVALS=TVALS+TVAL(I)
            TERSS=TERSS+TERS(I)
         ENDDO
         ROOT= DSQRT(2d0*TERSS)
         IF(ROOT .LE. ABSTOL .OR. ROOT .LE. RELTOL*DABS(TVALS)) GO TO 9
         IF(NTER .EQ. NDIM) GO TO 9
         BIGE=TERS(1)
         IBIG=1
         DO I = 2,NTER
            IF(TERS(I) .GT. BIGE) THEN
               BIGE=TERS(I)
               IBIG=I
            ENDIF
         ENDDO
         NTER=NTER+1
         XHI(NTER)=XHI(IBIG)
         XNEW=HF*(XLO(IBIG)+XHI(IBIG))
         XHI(IBIG)=XNEW
         XLO(NTER)=XNEW
         CALL DGS56P(F,XLO(IBIG),XHI(IBIG),TVAL(IBIG),TE)
         TERS(IBIG)=TE**2
         CALL DGS56P(F,XLO(NTER),XHI(NTER),TVAL(NTER),TE)
         TERS(NTER)=TE**2
      ENDDO
    9 RES=TVALS
      ERR=ROOT
      END SUBROUTINE DADAPT
!***************************************************************

!***************************************************************
      SUBROUTINE DGS56P(F,A,B,RES,ERR)
      IMPLICIT none

      integer:: I
      real(8),intent(in):: A,B
      real(8),intent(out):: RES,ERR
      real(8):: F, RANG, E5, E6
      external F
      real(8),parameter:: R1=1d0, HF=R1/2d0
 
      real(8),DIMENSION(5):: X5,W5
      real(8),DIMENSION(6):: X6,W6
  
      DATA (X5(I),W5(I),I=1,5)                            &
        / 4.6910077030668004D-02, 1.1846344252809454D-01, &
          2.3076534494715846D-01, 2.3931433524968324D-01, &
          5.0000000000000000D-01, 2.8444444444444444D-01, &
          7.6923465505284154D-01, 2.3931433524968324D-01, &
          9.5308992296933200D-01, 1.1846344252809454D-01  &
        /
 
      DATA (X6(I),W6(I),I=1,6)                            &
        / 3.3765242898423989D-02, 8.5662246189585178D-02, &
          1.6939530676686775D-01, 1.8038078652406930D-01, &
          3.8069040695840155D-01, 2.3395696728634552D-01, &
          6.1930959304159845D-01, 2.3395696728634552D-01, &
          8.3060469323313225D-01, 1.8038078652406930D-01, &
          9.6623475710157601D-01, 8.5662246189585178D-02  &
        /
      RANG=B-A
      E5=0d0
      E6=0d0
      do I = 1,5
         E5=E5+W5(I)*F(A+RANG*X5(I))
         E6=E6+W6(I)*F(A+RANG*X6(I))
      enddo
      E6=E6+W6(6)*F(A+RANG*X6(6))
      RES=HF*(E6+E5)*RANG
      ERR=DABS((E6-E5)*RANG)
      END SUBROUTINE DGS56P
!****************************************************************

!****************************************************************
      real(8) FUNCTION DGAUSS(F,A,B,EPS)
      implicit none

      real(8), intent(in)   :: A, B, EPS
      integer:: i
      real(8), DIMENSION(12):: W, X
      real(8), parameter    :: Z1=1d0, HF=Z1/2d0, CST=5D-3*Z1
      real(8)               :: F, CONST, C1, C2, S8, U, S16, AA, BB, H

      DATA X( 1) /9.6028985649753623D-1/, W( 1) /1.0122853629037626D-1/
      DATA X( 2) /7.9666647741362674D-1/, W( 2) /2.2238103445337447D-1/
      DATA X( 3) /5.2553240991632899D-1/, W( 3) /3.1370664587788729D-1/
      DATA X( 4) /1.8343464249564980D-1/, W( 4) /3.6268378337836198D-1/
      DATA X( 5) /9.8940093499164993D-1/, W( 5) /2.7152459411754095D-2/
      DATA X( 6) /9.4457502307323258D-1/, W( 6) /6.2253523938647893D-2/
      DATA X( 7) /8.6563120238783174D-1/, W( 7) /9.5158511682492785D-2/
      DATA X( 8) /7.5540440835500303D-1/, W( 8) /1.2462897125553387D-1/
      DATA X( 9) /6.1787624440264375D-1/, W( 9) /1.4959598881657673D-1/
      DATA X(10) /4.5801677765722739D-1/, W(10) /1.6915651939500254D-1/
      DATA X(11) /2.8160355077925891D-1/, W(11) /1.8260341504492359D-1/
      DATA X(12) /9.5012509837637440D-2/, W(12) /1.8945061045506850D-1/

      H=0d0
      IF(B .EQ. A) GO TO 99
      CONST=CST/DABS(B-A)
      BB=A
    1 AA=BB
      BB=B
    2 C1=HF*(BB+AA)
      C2=HF*(BB-AA)
      S8=0d0
      DO I = 1,4
         U=C2*X(I)
         S8=S8+W(I)*(F(C1+U)+F(C1-U))
      enddo
      S16=0d0
      DO I = 5,12
         U=C2*X(I)
         S16=S16+W(I)*(F(C1+U)+F(C1-U))
      enddo
      S16=C2*S16
      IF(DABS(S16-C2*S8) .LE. EPS*(1+DABS(S16))) THEN
         H=H+S16
         IF(BB .NE. B) GO TO 1
      ELSE
         BB=C1
         IF(1+CONST*DABS(C2) .NE. 1) GO TO 2
         H=0d0
         write(*,'(''DGAUSS...TOO HIGH ACCURACY REQUIRED'')')
         GO TO 99
      END IF

   99 DGAUSS=H
      END FUNCTION DGAUSS
!***************************************************************

!***************************************************************
! $ H = \int_A^B F(x) dx $
  Subroutine GAUSS(F,A,B,EPS,H)
     implicit none      
     integer:: i
     real(8), intent(in)   :: A, B, EPS
     real(8), intent(out)  :: H
     real(8), DIMENSION(12):: W, X
     real(8), parameter    :: Z1=1d0, HF=Z1/2d0, CST=5D-3*Z1
     real(8)               :: F, CONST, C1, C2, S8, U, S16, AA, BB

      DATA X( 1) /9.6028985649753623D-1/, W( 1) /1.0122853629037626D-1/
      DATA X( 2) /7.9666647741362674D-1/, W( 2) /2.2238103445337447D-1/
      DATA X( 3) /5.2553240991632899D-1/, W( 3) /3.1370664587788729D-1/
      DATA X( 4) /1.8343464249564980D-1/, W( 4) /3.6268378337836198D-1/
      DATA X( 5) /9.8940093499164993D-1/, W( 5) /2.7152459411754095D-2/
      DATA X( 6) /9.4457502307323258D-1/, W( 6) /6.2253523938647893D-2/
      DATA X( 7) /8.6563120238783174D-1/, W( 7) /9.5158511682492785D-2/
      DATA X( 8) /7.5540440835500303D-1/, W( 8) /1.2462897125553387D-1/
      DATA X( 9) /6.1787624440264375D-1/, W( 9) /1.4959598881657673D-1/
      DATA X(10) /4.5801677765722739D-1/, W(10) /1.6915651939500254D-1/
      DATA X(11) /2.8160355077925891D-1/, W(11) /1.8260341504492359D-1/
      DATA X(12) /9.5012509837637440D-2/, W(12) /1.8945061045506850D-1/

      H=0d0
      IF(B .EQ. A) return
      CONST=CST/DABS(B-A)
      BB=A
    1 AA=BB
      BB=B
    2 C1=HF*(BB+AA)
      C2=HF*(BB-AA)
      S8=0d0

      DO I = 1,4
         U=C2*X(I)
         S8=S8+W(I)*(F(C1+U)+F(C1-U))
      enddo
      S16=0d0

      DO I = 5,12
         U=C2*X(I)
         S16=S16+W(I)*(F(C1+U)+F(C1-U))
      enddo
      S16=C2*S16

      IF(DABS(S16-C2*S8) .LE. EPS*(1d0 + DABS(S16))) THEN
         H = H+S16
         IF(BB .NE. B) GO TO 1
      ELSE
         BB = C1
         IF(1d0 + CONST*DABS(C2) .NE. 1d0) GO TO 2
         H=0d0
         write(*,'(''GAUSS...TOO HIGH ACCURACY REQUIRED'')')
         return
      END IF

  END SUBROUTINE GAUSS
!***************************************************************

!***********************************************************************
  SUBROUTINE QUANC8(FUN,A,B,ABSERR,RELERR,RESULT,ERREST,NOFUN,FLAG)
!
      REAL(8) :: FUN
      REAL(8), intent(in) :: A,B,ABSERR,RELERR
      REAL(8), intent(out) :: RESULT, ERREST,FLAG
      INTEGER, intent(out) :: NOFUN
!
      REAL(8) :: W0, W1, W2, W3, W4, STEP   &
               , AREA, X0, F0, STONE, COR11 &
               , TEMP, QPREV, QNOW, QDIFF   &
               , QLEFT, ESTERR, TOLERR
      REAL(8), DIMENSION(31)   :: QRIGHT
      REAL(8), DIMENSION(16)   :: F, X
      REAL(8), DIMENSION(8,30) ::FSAVE, XSAVE
      INTEGER :: LEVMIN,LEVMAX,LEVOUT,NOMAX,NOFIN,LEV,NIM,I,J
!
      LEVMIN = 1
      LEVMAX = 30
      LEVOUT = 6
      NOMAX  = 5000
      NOFI   = NOMAX-8*(LEVMAX-LEVOUT+2**(LEVOUT+1))
!
      W0 =  3956d0/14175d0
      W1 =  23552d0/14175d0
      W2 = -3712d0/14175d0
      W3 =  41984d0/14175d0
      W4 = -18160d0/14175d0
!
      FLAG   = 0d0
      RESULT = 0d0
      COR11  = 0d0
      ERREST = 0d0
      AREA   = 0d0
      NOFUN  = 0
      IF(A.EQ.B) RETURN
!
      LEV   = 0
      NIM   = 1
      X0    = A
      X(16) = B
      QPREV = 0d0
      F0    = FUN(X0)
      STONE = (B-A)/16d0
      X(8)  = (X0+X(16))/2d0
      X(4)  = (X0+X(8))/2d0
      X(12) = (X(8)+X(16))/2d0
      X(2)  = (X0+X(4))/2d0
      X(6)  = (X(4)+X(8))/2d0
      X(10) = (X(8)+X(12))/2d0
      X(14) = (X(12)+X(16))/2d0
      DO J = 2,16,2
         F(J)=FUN(X(J))
      ENDDO
      NOFUN = 9
!
30    X(1) = (X0 + X(2))/2d0
      F(1) = FUN(X(1))
      DO J = 3,15,2
         X(J) = (X(J-1) + X(J+1))/2d0
2        F(J) = FUN(X(J))
      ENDDO
      NOFUN = NOFUN + 8
      STEP  = (X(16)-X0)/16d0
      QLEFT = (W0*(F0 + F(8)) + W1*(F(1) + F(7)) + W2*(F(2) + F(6)) &
             + W3*(F(3) + F(5)) + W4*F(4))*STEP
      QRIGHT(LEV+1) = (W0*(F(8) + F(16)) + W1*(F(9) + F(15))  &
                    + W2*(F(10) + F(14)) + W3*(F(11) + F(13)) &
                    + W4*F(12))*STEP
      QNOW  = QLEFT + QRIGHT(LEV+1)
      QDIFF = QNOW  - QPREV
      AREA  = AREA  + QDIFF
!
      ESTERR = DABS(QDIFF)/1023d0
      TOLERR = DMAX1(ABSERR,RELERR*ABS(AREA))*(STEP/STONE)
      IF(LEV.LT.LEVMIN) GO TO 50
      IF(LEV.GE.LEVMAX) GO TO 62
      IF(NOFUN.GT.NOFIN) GO TO 60
      IF(ESTERR.LE.TOLERR) GO TO 70
!
50    NIM=2*NIM
      LEV=LEV+1
!
      DO I=1,8
         FSAVE(I,LEV)=F(I+8)
         XSAVE(I,LEV)=X(I+8)
      ENDDO
!
      QPREV = QLEFT
      DO I =  1,8
         J = -I
         F(2*J+18) = F(J+9)
         X(2*J+18) = X(J+9)
      ENDDO
      GO TO 30
!
60    NOFIN  = 2*NOFIN
      LEVMAX = LEVOUT
      FLAG=FLAG+(B-X0)/(B-A)
      GO TO 70
!
62    FLAG = FLAG + 1d0
!
70    RESULT = RESULT + QNOW
      ERREST = ERREST + ESTERR
      COR11  = COR11  + QDIFF/1023d0
!
72    IF(NIM.EQ.2*(NIM/2)) GO TO 75
      NIM = NIM/2
      LEV = LEV - 1
      GO TO 72
75    NIM = NIM+1
      IF(LEV.LE.0) GO TO 80
!
      QPREV = QRIGHT(LEV)
      X0 = X(16)
      F0 = F(16)
      DO I=1,8
         F(2*I)=FSAVE(I,LEV)
         X(2*I)=XSAVE(I,LEV)
      ENDDO
      GO TO 30
!
80    RESULT = RESULT + COR11
!
      IF(ERREST.EQ.0d0) RETURN
82    TEMP = DABS(RESULT)+ERREST
      IF(TEMP.NE.DABS(RESULT)) RETURN
      ERREST = 2d0*ERREST
      GO TO 82
 END SUBROUTINE QUANC8
!*******************************************************************

!***************************************************************!
!                      INTERPOLATION                            ! 
!***************************************************************!

!*******************************************************************
! Linear interpolation and extrapolation of function F(x);
! ya=F(xa); yb=F(xb) 
  real(8) Function FLin(x,xa,xb,ya,yb)
     implicit none
     real(8), intent(in) :: x, xa, xb, ya, yb
     if(xa.eq.xb) then
        print*, 'in FLin xa=xb ! Stop!'
        stop
     endif 
     Flin = ya + (yb-ya)/(xb-xa)*(x-xa)  
  end Function FLin
!*******************************************************************


!*******************************************************************
!
  real(8) FUNCTION DIVDIF(F,X,N,Z,M)
     implicit none
     integer,intent(in)             :: N
     real(8),intent(in),dimension(N):: F,X
     real(8),intent(in)             :: Z 
     real(8),dimension(11)          :: COF, ARG
     integer                        :: MMAX = 10
     real(8)                        :: ZERO = 0d0
     real(8)                        :: COF1, COF2, SUM
     integer                        :: M , MM, IPOINT, I , II, J   &
                                     , IL, IU, JL    , JU, I1, I2  &
                                     , INDEX, JNDEX
!
!     CERN LIBRARY PROGRAM NO
!
!     BASIC VERSION APRIL 6971.
!
!     PURPOSE = TO INTERPOLATE IN TABLE OF GIVEN FUNCTION VALUES WHICH
!               ARE STORED AFTER INCREASING OR DECREASING VALUES OF THE
!               ARGUMENTS.
!
!     PARAMETERS
!
!     F       = THE ARRAY OF THE GIVEN FUNCTION VALUES.F(K)=F(X(,)).
!     X       = THE ARRAY OF GIVEN ARGUMENTS.
!     N       = DIMENSION OF THE ARRAYS F AND X,I.E.THE NUMBER OF POINTS
!               TABLE.
!     Z       = ARGUMENT FOR WHICH THE INTERPOLATION IS WANTED.
!     M       = ORDER OF INTERPOLATION.

!
!     INTERNAL PARAMETERS.
!
!     MMAX    = THE MAXIMUM ORDER OF INTERPOLATION PERMITTED.
!     COF     = AN INTERNAL ARRAY DIMENSIONED COF(MMAX+1) USED FOR STO-
!               RING THE COMPUTED COEFFICIENTS IN THE INTERPOLATION FOR-
!               MULA.
!     ARG     = AN INTERNAL ARRAY DIMENSIONED ARG(MMAX+1) USED FOR STO-
!               RING THE ABSCISSAE VALUES USED IN THE INTERPOLATION.
!
 1000 IF ((Z-X(1))*(X(N)-Z).GE.ZERO) GO TO 1030
!
!     Z-VALUE OUTSIDE RANGE,PRINT ERROR MESSAGE.
!
 1020 PRINT 10 , Z
      stop
      DIVDIF = ZERO
      RETURN
!
 1030 IF ((M.LE.(N-1)).AND.(M.LE.MMAX)) GO TO 1040
      MM=M
      IF (M.GT.(N-1)) M=N-1
      IF (M.GT.MMAX) M=MMAX
!
!     REQUIRED ORDER OF INTERPOLATION TOO HIGH.PRINT ERROR MESSAGE AND
!     REDUCE ORDER.
!
      PRINT 20 , MM , M
!
!     START ACTUAL CALCULATION.
!
!     COMPUTE POINTER,IPOINT,FOR THE LEFT BOUNDARY OF THE INTERVAL IN
!     WHICH WE HAVE Z I.E. Z IN THE INTERVAL X(IPOINT),X(IPOINT+1).
!
 1040 COF1 = Z - X(1)
      DO I = 2,N
         IPOINT = I-1
         COF2 = Z - X(I)
         IF (COF1*COF2.LE.ZERO) GO TO 1060
         COF1 = COF2
      ENDDO 
!
!     CONSTRUCT TABLE TO BE USED IN THE INTERPOLATION.
!
 1060 IL = IPOINT
      IU = IL+1
      JL = 1
      JU = 1
      MM = M+1
      DO I=1,MM
         I1=1
         I2=1
         IF ((JL.EQ.0).OR.(JU.EQ.0)) GO TO 1070
         COF1=DABS(Z-X(IL))
         COF2=DABS(X(IU)-Z)
         IF (COF1.GT.COF2) I1=0
         IF (I1.EQ.1) I2=0
 1070    IF ((JL.EQ.0).OR.(I1.EQ.0)) II=IU
         IF ((JU.EQ.0).OR.(I2.EQ.0)) II=IL
         ARG(I) = X(II)
         COF(I) = F(II)
         IF ((JL.EQ.1).AND.(I1.EQ.1)) IL=IL-1
         IF ((JU.EQ.1).AND.(I2.EQ.1)) IU=IU+1
         IF (IL.LT.1) JL=0
         IF (IU.GT.N) JU=0
      ENDDO
!
      DO I = 1,M
         DO J = I,M
            INDEX = M + 1 + I - J
            JNDEX = INDEX - I
            COF(INDEX) = (COF(INDEX)-COF(INDEX-1)) &
          / (ARG(INDEX)-ARG(JNDEX))
         ENDDO
      ENDDO
!
      SUM = COF(M+1)
      DO I = 1,M
         INDEX=M+1-I
         SUM = (Z-ARG(INDEX))*SUM+COF(INDEX)
      ENDDO
!
      DIVDIF = SUM
!
 10   FORMAT(//,5X,'*** ERROR MESSAGE FUNCTION DIVDIF *** , ARGUMENT Z = ' &
          , D21.14,' OUTSIDE RANGE.',//)
 20   FORMAT(//,5X,'*** ERROR MESSAGE FUNCTION DIVDIF *** , ORDER OF INTERPOLATION M =' &
          , I3,' TOO HIGH,M IS REDUCED TO' ,I2,//)
!
  END FUNCTION DIVDIF
!*******************************************************************

!*******************************************************************
!  2-dim linear interpolaton
!  simultaneously linearly interpolates L arrays of dim. NxM
!  result - vector of interpolated values of dim. L
  SUBROUTINE EXTVEC(X,Y,F,N,M,L,XX,YY,FF)
     implicit none
     integer,intent( in)                  :: N, M, L 
     real(8),intent( in)                  :: XX, YY   
     real(8),intent( in), dimension(N)    :: X 
     real(8),intent( in), dimension(M)    :: Y 
     real(8),intent( in), dimension(N,M,L):: F 
     real(8),intent(out), dimension(L)    :: FF 
     integer                              :: I, J, II, JJ, K1
     real(8)                              :: R1, R2, F1, F2 

!-----------------------------
     if(XX.le.X(1)) then
        I=1
     elseif(XX.ge.X(N)) then
        I=N-1
     else  
        DO II=1,N
           IF (XX.lt.X(II)) exit
        ENDDO
        I=II-1
     endif
!-----------------------------
     if(YY.le.Y(1)) then
        J=1
     elseif(YY.ge.Y(M)) then
        J = M-1
     else  
        DO JJ=1,M
           IF (YY.lt.Y(JJ)) exit
        ENDDO
        J = JJ-1
     endif
!-----------------------------
     R1 = X(I+1) - X(I)
     IF (R1.eq.0d0) then
        write(*,'(" ERROR IN DATA",10X,"I=",I4,5X,"X(I)=",d10.3, &
           5X,"X(I+1)=",d10.3)') I, X(I), X(I+1)
        stop
     ENDIF

     R2 = Y(J+1) - Y(J)
     IF (R2.eq.0d0) then
        write(*,'(" ERROR IN DATA",10X,"J=",I4,5X,"Y(J)=",d10.3, &
           5X,"Y(J+1)=",d10.3)') J, Y(J), Y(J+1)
        stop
     ENDIF

     R1 = (XX-X(I))/R1

     DO K1 = 1,L
        F1 = F(I,J+1,K1) + (F(I+1,J+1,K1) - F(I,J+1,K1))*R1
        F2 = F(I,J,K1) + R1*(F(I+1,J,K1) - F(I,J,K1))
        FF(K1) = F2 + (YY - Y(J))*(F1 - F2)/R2
     ENDDO

  END SUBROUTINE EXTVEC
!*******************************************************************

!*******************************************************************
! Find lower position 'i' in the 1-dim array XX(N)
!         XX(i) <= x < XX(i+1)
! May be used for following interpolation
  Subroutine FindPos(N, XX, x, i)
     implicit none
     integer, intent(in) :: N
     real(8), intent(in) :: x
     real(8), intent(in), dimension(N) :: XX
     integer, intent(out) :: i 

     integer :: itemp, ilow, iup

     i = 1
     if((x-XX(1))*(x-XX(N)).gt.0d0) then
        print*, ' Attention, x is out of range in FindPos!!!'
        stop
     endif
    
	 if(x.eq.xx(N)) then
        i = n
     else
	    if(N.lt.3) then
           i = 1
           return
         endif

        ilow = 1
        iup  = N
!------- 
        if(x.lt.XX(2)) then
           i=1
           return
        endif
1       itemp = (ilow + iup)/2
        if(x.ge.XX(itemp).and.x.lt.XX(itemp+1)) then
           i = itemp
           return
        endif
        if(x.lt.XX(itemp+1)) then   
           iup = itemp+1
        else 
           ilow = itemp
        endif
        goto 1
     endif
  end Subroutine FindPos
!*******************************************************************

!***************************************************************!
!             ORDINARY DIFFERENTIAL EQUATION                    !
!***************************************************************!
!     Based on a modification of the Runge-Kutta method suggested
!     by Merson. See G.N. Lance, Numerical Methods for High speed
!     Computers, Iliffe & Sons, London 1960, pp. 56-57

      SUBROUTINE DDEQMR(N,XA,XZ,Y,H0,EPS,SUB,W)
      implicit none 
  
      real(8),intent(in)    :: XA, XZ, EPS
      integer,intent(in)    :: N
      real(8),dimension(*)  :: Y
      real(8),dimension(N,*):: W    

      CHARACTER(80):: ERRTXT
      LOGICAL:: LER, LFN

      real(8):: delta = 1d-14, H0, X
      real(8),parameter:: Z1 = 1d0, R2 = Z1/2d0, R3 = Z1/3d0 &
                        , R4 = 3d0*Z1/8d0, R5 = 3d0*Z1/2d0   &
                        , R6 = 9d0*Z1/2d0, R7 = 4d0*Z1/3d0, R0 = Z1/32d0 
      real(8):: DELTAX, EPS5, EPS0, H1, SGH, HH, S2, S3, S7, X1, X2, X3
      integer:: I

      IF(N .LT. 1 .OR. XA .EQ. XZ .OR. H0 .EQ. 0D0) RETURN
      DELTAX=DELTA*DABS(XZ-XA)
      EPS5=5d0*DABS(EPS)
      EPS0=R0*EPS5
      X=XA
      H1=DSIGN(DABS(H0),XZ-XA)
      SGH=DSIGN(Z1,H1)

   12 IF(SGH*(X+H1-XZ) .LT. 0) THEN
       HH=H1
       H0=H1
       LFN=.FALSE.
      ELSE
       HH=XZ-X
       IF(DABS(HH) .LT. DELTAX) THEN
        DO 10 I = 1,N
   10   Y(I)=W(I,6)
        RETURN
       END IF
       LFN=.TRUE.
      END IF
      S2=R2*HH
      S3=R3*HH
      S7=R7*HH
      X1=X+HH
      X2=X+S2
      X3=X+S3

      CALL SUB(X,Y,W(1,1))
      DO 1 I = 1,N
      W(I,1)=S3*W(I,1)
    1 W(I,6)=Y(I)+W(I,1)

      CALL SUB(X3,W(1,6),W(1,2))
      DO 2 I = 1,N
      W(I,2)=S3*W(I,2)
    2 W(I,6)=Y(I)+R2*(W(I,1)+W(I,2))

      CALL SUB(X3,W(1,6),W(1,3))
      DO 3 I = 1,N
      W(I,3)=S3*W(I,3)
      W(I,2)=3*W(I,3)
    3 W(I,6)=Y(I)+R4*(W(I,1)+W(I,2))

      CALL SUB(X2,W(1,6),W(1,4))
      DO 4 I = 1,N
      W(I,4)=S7*W(I,4)
    4 W(I,6)=Y(I)+R5*(W(I,1)-W(I,2)+W(I,4))

      CALL SUB(X1,W(1,6),W(1,5))
      DO 5 I = 1,N
      W(I,5)=S3*W(I,5)
    5 W(I,6)=Y(I)+R2*(W(I,1)+W(I,4)+W(I,5))

      DO 8 I = 1,N
      W(I,2)=DABS(W(I,1)-R6*W(I,3)+W(I,4)-R2*W(I,5))
      W(I,1)=DABS(W(I,6))
      IF(W(I,2) .GT. EPS5*W(I,1)) THEN
       H1=R2*HH
       IF(DABS(H1) .LT. DELTAX) THEN
        WRITE(ERRTXT,101) X
        RETURN
       END IF
       GO TO 12
      END IF
    8 CONTINUE
      LER=.TRUE.
      DO 7 I = 1,N
    7 LER=LER .AND. W(I,2) .LT. EPS0*W(I,1)
      DO 9 I = 1,N
    9 Y(I)=W(I,6)
      IF(LER) THEN
       H0=H1+H1
       H1=HH+HH
      END IF
      IF(LFN) RETURN
      X=X1
      GO TO 12
  101 FORMAT('TOO HIGH ACCURACY REQUIRED NEAR  X = ',1P,D15.8)
      END SUBROUTINE DDEQMR
!***************************************************************

!*******************************************************************
   real(8) Function ERF(X)
!
!=========================================
!  Purpose: Compute error function erf(x)
!  Input:   x   --- Argument of erf(x)
!  Output:  ERR --- erf(x)
!========================================
!  
       IMPLICIT NONE
       real(8),intent(in)::  X
       real(8),parameter  :: EPS = 1.0D-15 &
                           , PI  =3.141592653589793D0
       real(8):: X2, ER, R, C0
       integer:: K

        X2=X*X

        IF (DABS(X).LT.3.5D0) THEN
           ER = 1.0D0
           R  = 1.0D0
           DO K=1,50
              R=R*X2/(K+0.5D0)
              ER=ER+R
              IF (DABS(R).LE.DABS(ER)*EPS) GO TO 15
           ENDDO
15         C0=2.0D0/DSQRT(PI)*X*DEXP(-X2)
           ERF = C0*ER
        ELSE
           ER=1.0D0
           R=1.0D0
           DO K=1,12
              R=-R*(K-0.5D0)/X2
              ER=ER+R
           ENDDO 
           C0=DEXP(-X2)/(DABS(X)*DSQRT(PI))
           ERF = 1.0D0-C0*ER
           IF (X.LT.0d0) ERF = -ERF
        ENDIF
        END FUNCTION ERF
!*******************************************************************

! ----------------------------------------------------------------------
   real(8) FUNCTION BESSK0(X)
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
!     POUR TOUT X REEL NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
! ----------------------------------------------------------------------
      implicit none
      real(8),intent(in):: x
      REAL(8):: Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7

      DATA P1,P2,P3,P4,P5,P6,P7                        &
      / -0.57721566D0,0.42278420D0,0.23069756D0,       &
         0.3488590D-1,0.262698D-2,0.10750D-3,0.74D-5/
      
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7                        &
      / 1.25331414D0,-0.7832358D-1,0.2189568D-1,       & 
        -0.1062446D-1,0.587872D-2,-0.251540D-2,0.53208D-3/
      
      IF(X.EQ.0.D0) THEN
         BESSK0=1.D30
         RETURN
      ENDIF

      IF(X.LE.2.D0) THEN
         Y      = X*X/4.D0
         AX     = -DLOG(X/2.D0)*BESSI0(X)
         BESSK0 = AX+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         Y  = (2.D0/X)
         AX = DEXP(-X)/DSQRT(X)
         BESSK0 = AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF

      END FUNCTION BESSK0
! ----------------------------------------------------------------------
   real(8) FUNCTION BESSK1(X)
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
!     POUR TOUT X REEL POSITF NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
! ----------------------------------------------------------------------
      implicit none
      real(8),intent(in):: x
      REAL(8):: Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7
      
      DATA P1,P2,P3,P4,P5,P6,P7           &
      / 1.D0,0.15443144D0,-0.67278579D0,  &
       -0.18156897D0,-0.1919402D-1,-0.110404D-2,-0.4686D-4/
      
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7                  &
      / 1.25331414D0,0.23498619D0,-0.3655620D-1, &
        0.1504268D-1,-0.780353D-2,0.325614D-2,-0.68245D-3/

      IF(X.EQ.0.D0) THEN
        BESSK1 = 1.D32
        RETURN
      ENDIF
      
      IF(X.LE.2.D0) THEN
         Y  = X*X/4.D0
         AX = DLOG(X/2.D0)*BESSI1(X)
         BESSK1=AX+(1.D0/X)*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         Y  = (2.D0/X)
         AX = DEXP(-X)/DSQRT(X)
         BESSK1 = AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      ENDIF

   END FUNCTION BESSK1
! ----------------------------------------------------------------------
   real(8) FUNCTION BESSK2(X)
      implicit none
      real(8),intent(in):: x

      BESSK2 = BESSK0(x) + 2d0*BESSK1(x)/x

   end FUNCTION BESSK2
! ----------------------------------------------------------------------
! Auxiliary Bessel functions for N=0, N=1
   real(8) FUNCTION BESSI0(X)
      implicit none
      real(8),intent(in):: x
      REAL(8):: Y,P1,P2,P3,P4,P5,P6,P7,            &
                Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX

      DATA P1,P2,P3,P4,P5,P6,P7                    &
      / 1.D0,3.5156229D0,3.0899424D0,1.2067429D0,  &
        0.2659732D0,0.360768D-1,0.45813D-2/

      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9              &
      / 0.39894228D0,0.1328592D-1, .225319D-2,     &
       -0.157565D-2,0.916281D-2,-0.2057706D-1,     &
        0.2635537D-1,-0.1647633D-1,0.392377D-2/

      IF(ABS(X).LT.3.75D0) THEN
         Y = (X/3.75D0)**2
         BESSI0 = P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
         AX = DABS(X)
         Y  = 3.75D0/AX
         BX = DEXP(AX)/DSQRT(AX)
         AX = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
         BESSI0 = AX*BX
      ENDIF

   END FUNCTION BESSI0
! ----------------------------------------------------------------------
   real(8) FUNCTION BESSI1(X)
      implicit none
      real(8),intent(in):: x
      REAL(8):: Y,P1,P2,P3,P4,P5,P6,P7,                     &
                Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX

      DATA P1,P2,P3,P4,P5,P6,P7                             &
     / 0.5D0,0.87890594D0,0.51498869D0,                    &
        0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
      
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9                       &
      / 0.39894228D0,-0.3988024D-1, &
       -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, &
       -0.2895312D-1,0.1787654D-1,-0.420059D-2/
      
      IF(ABS(X).LT.3.75D0) THEN
         Y      = (X/3.75D0)**2
         BESSI1 = X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         AX = DABS(X)
         Y  = 3.75D0/AX
         BX = DEXP(AX)/DSQRT(AX)
         AX = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))
         BESSI1 = AX*BX
      ENDIF

   ENDFUNCTION BESSI1
! ----------------------------------------------------------------------

!*******************************************************************
!---------------------------------!
! Try Integrator of Vadim Naumov  !
!---------------------------------!
!   Program Naum
!      use gazlib,only: dgauss, mulset, mulint, multot, DADAPT 
!      implicit none
!
!      external FNaum, FNaum1
!
!      real(8):: xl, xu, relerr, res, Pi, FNaum, FNaum1
!      real(8):: resmul, Aerr, Err
!
!      integer,parameter:: MinCal = 100  &
!              , NN1 = 1, NN2 = 2
!
!      real(8),dimension(2):: Xlow, Xupp
!
!      COMMON /MulLim/ Xlow, Xupp
!=====================================================
!
!      Pi = dacos(-1d0)
!      xl = 0d0  ! -Pi/2d0
!      xu = 1d0  !  Pi/2d0
!      relerr = 1d-14
!      Aerr = 1d-15
!
!      resmul = 0d0
!
!-----------------------------------------------------------
!      call DADAPT(FNaum1,xl,xu,1,relerr,Aerr,RES,ERR)
!      write(*,'(A,E20.13)') '  Using DADAPT, res =', res
!
!      res = dgauss(FNaum1,xl,xu,relerr)
!      write(*,'(A,E20.13)') '  Using dgauss, res =', res
!-----------------------------------------------------------
!      Xlow(1) = xl
!      Xupp(1) = xu
!      Xlow(2) = xl
!      Xupp(2) = xu
!
!      call mulset(FNaum1,Res,RelErr,Aerr,MinCal,NN1)
!      call mulint(FNaum1,Res)
!      write(*,'(A,E20.13)') '  Using mulN1 , res =', res
!
!     call mulset(FNaum,Resmul,RelErr,Aerr,MinCal,NN2)
!     call multot(FNaum,Resmul,RelErr,Aerr,MinCal)
!      write(*,'(A,E20.13)') '  Using mulN2 , res =', Resmul
!
!      stop
!
!   end Program Naum
!
!---------------------------------
!   real(8) Function FNaum(x)
!      real(8),dimension(2) :: x
!      FNaum = dexp(-1d0/x(1))
!   end Function FNaum 
!
!---------------------------------
!   real(8) Function FNaum1(x)
!      real(8):: x
!      FNaum1 = dexp(-1d0/x)
!   end Function FNaum1
!---------------------------------


!*******************************************************************
! Updated from Vadim Naumov 
!***********************************************************************
! ******************************************************************** *
! *                                                                  * *
! *   MM      MM     AA     TTTTTTTT  HH    HH             QQQQQQ    * *
! *   MMM    MMM    A  A       TT     HH    HH            QQ    QQ   * *
! *   MM M  M MM   A    A      TT     HH    HH            QQ    QQ   * *
! *   MM  MM  MM  AAAAAAAA     TT     HHHHHHHH            QQ    QQ   * *
! *   MM      MM  AA    AA     TT     HH    HH            QQ    QQ   * *
! *   MM      MM  AA    AA     TT     HH    HH            QQ  Q QQ   * *
! *   MM      MM  AA    AA     TT     HH    HH  =========  QQQQQQ    * *
! *                                                             QQ   * *
! *         General mathematical procedures: QUADRATURES             * *
! *                                                                  * *
! ******************************************************************** *
!                                                                      *
!  ##################################################################  *
!  #                                                                #  *
!  #     Vadim Naumov, INFN, Sezione di Firenze and LTP, ISU        #  *
!  #                                                                #  *
!  # Version of August 15, 1997  for MS 32-bit FORTRAN PowerStation #  *
!  #                                                                #  *
!  ##################################################################  *
!                                                                      *
! ******************************************************************** *
! *                                                                  * *
! *   MM      MM            ll          -------------------------    * *
! *   MMM    MMM            ll         | ADAPTIVE QUADRATURE FOR |   * *
! *   MM M  M MM  uu    uu  ll         | MULTIPLE INTEGRALS OVER |   * *
! *   MM  MM  MM  uu    uu  ll         |   RECTANGULAR REGIONS   |   * *
! *   MM      MM  uu    uu  ll         | by Genz-Malik algorithm |   * *
! *   MM      MM  uu    uu  ll         |    (revized version)    |   * *
! *   MM      MM    uuuuuu   llllll     -------------------------    * *
! *                                                                  * *
! *                                                                  * *
! *   ------------------------------------------------------------   * *
! *                                                                  * *
! *                       T H E    S O U R C E                       * *
! *                                                                  * *
! *    --------                                            ------    * *
! *   | RADMUL |          CERN Program Library            | D120 |   * *
! *    --------                                            ------    * *
! *                                                                  * *
! *   Authors: A.C.Genz & A.A.Malik          Library: MATHLIB        * *
! *   Submitter: K.S.Kolbig                  Submitted: 15.11.1995   * *
! *                                                                  * *
! *   References:                                                    * *
! *                                                                  * *
! *   [1] A.C.Genz & A.A.Malik, Remarks on algorithm 006: An adap-   * *
! *       tive algorithm for numerical integration over N-dimensi-   * *
! *       onal rectangular region, J. Comput. Appl. Math. 6 (1980)   * *
! *       295-302.                                                   * *
! *   [2] A. van Doren & L. de Ridder,  An adaptive algorithm  for   * *
! *       numerical integration over N-dimensional cube, J.Comput.   * *
! *       Appl. Math. 2 (1976) 207-217.                              * *
! *                                                                  * *
! *   ------------------------------------------------------------   * *
! *   E N T R I E S:                                                 * *
! *   MulSet(Fun,Res,RelErr,Aerr,MinCal,IntDim)                      * *
! *   MulTot(Fun,Res,RelErr,Aerr,MinCal)                             * *
! *   MulInt(Fun,Res)                                                * *
! *   MulInf                                                         * *
! *   ------------------------------------------------------------   * *
! *                                                                  * *
! ******************************************************************** *
!***********************************************************************
      SUBROUTINE MulSet(Fun,Res,RelErr,Aerr,MinCal,IntDim)
!***********************************************************************
!
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      SAVE
      LOGICAL*2 Flag

!     ------------------------------------------------------------------

!     THE FOLLOWING THREE PARAMETERS MUST BE SPECIFIED BEFORE USING:

      PARAMETER(Length = 90000000)      ! > (2*MaxDim+3)(1+k)/2, k > 0
      PARAMETER(MaxDim =        2)      ! > 1 [MAX(2,IntDim)]
      PARAMETER(Nlog  =         1)      ! LogFile number (open in MAIN)

!     ------------------------------------------------------------------

      PARAMETER(Zero = 0, One = 1, Half = One/2)
      PARAMETER(V2 = 245.0d0/486,   V4 =   25.0d0/729  )
      PARAMETER(W2 = 980.0d0/6561,  W4 =  200.0d0/19683)
      PARAMETER(C0 =   4.0d0/19683, C1 =    1.0d0/729  )
      PARAMETER(C3 =   5.0d0/1458,  C5 = 6859.0d0/19683)
      PARAMETER(Y2=0.3585685828003181d00)  ! 3/SQRT(70)
      PARAMETER(Y4=0.9486832980505138d00)  ! 3/SQRT(10)
      PARAMETER(Y5=0.6882472016116853d00)  ! 3/SQRT(19)
!**** DIMENSION STORE(Length)              ! storage array
      ALLOCATABLE STORE(:)                 ! storage array
      DIMENSION CTR(MaxDim),WTH(MaxDim),VTH(MaxDim),X(MaxDim)

!     ------------------------------------------------------------------
      COMMON /MulLim/Xlow(MaxDim),Xupp(MaxDim)
!     the arrays are lower limits/upper limits of integration
!     ------------------------------------------------------------------

!**** INCLUDE 'Heap.fi'
!**** EQUIVALENCE (STORE(1),Heap(1))       ! memory saving
      EQUIVALENCE (CTR1,CTR(1)),(WTH1,WTH(1)),(WTH2,WTH(2))

!     ------------------------------------------------------------------
!        THIS FRAGMENT MAY BE REMOVED AS A UNIT WHEN DEBUGGED          !
         IF(IntDim.LT.1.OR.IntDim.GT.MaxDim.OR.MinCal.LT.1) THEN       !
         WRITE(Nlog,*) ' '                                             !
         WRITE(Nlog,*) '   ------------------------ '                  !
         WRITE(Nlog,*) '  | MISTAKE AT CALL MULSET |'                  !
         WRITE(Nlog,*) '   ------------------------ '                  !
         WRITE(Nlog,*) ' '                                             !
         WRITE(Nlog,*) '   MaxDim =',MaxDim                            !
         WRITE(Nlog,*) '   IntDim =',IntDim                            !
         WRITE(Nlog,*) '   MinCal =',MinCal                            !
         WRITE(Nlog,*) ' '                                             !
         STOP          ' MISTAKE AT CALL MULSET     '                  !
      endIF                                                            !
!     ------------------------------------------------------------------

      EPS = RelErr              ! "relative integration accuracy"
      Aerrr = Aerr

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!                                                                      c
!  COMMENTS # The authors of the method credited the parameter EPS as  c
!             being the "specified relative accuracy of integration".  c
!             Unfortunately, this is not the case,  especially in the  c
!             original REAL*4 version for VAX. Generally this parame-  c
!             ter (together with parameters Length and MinCal) provi-  c
!             des a control upon the accuracy and a means for minimi-  c
!             zation of CPU time and memory.  A good rule of thumb is  c
!             "THE REAL INTEGRATION ACCURACY RISES AS EPS DECREASES".  c
!             Usually, for "not-too-terrible" integrands the real ac-  c
!             curacy is BETTER than EPS but from a certain (dependent  c
!             of the integral) value,  further decrease of EPS leaves  c
!             the result invariant. It is therefore strongly recomme-  c
!             nded that  FOR EVERY NEW INTEGRAND AND/OR EVERY NEW SET  c
!             OF THE LIMITS OF INTEGRATION, THE INTEGRATION PROCEDURE  c
!             BE REPEATED WITH DIFFERENT EPS.                          c
!                                                                      c
!           # Use entry 'MulInt' for any series of integrals with the  c
!             common parameters RelErr, MinCal and IntDim.             c
!                                                                      c
!           # Use entry 'MulTot' if, with the same IntDim, the RelErr  c
!             and/or MinCal must be changed.                           c
!                                                                      c
!           # When the parameters  Length  [the length of array STORE  c
!             of working storage],  MinCal [the minimum number of in-  c
!             tegrand Fun calls requested to be allowed] and RelErr =  c
!             EPS are fitted closely,  the VAX REAL*16 version of the  c
!             routine gives very good results  (accuracy-time-memory)  c
!             for integrands with variation up to (and including) 150  c
!             ORDERS OF MAGNITUDE!                                     c
!                                                                      c
!           # Maximum multiplicity is formally unbounded but actually  c
!             it is limited by your computer. For instance, ALPHA-VAX  c
!             usually has no troubles with 20-fold integrals.          c
!                                                                      c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!     IF(Res==Zero) Res=Fun(X)       ! An optional possibility.

      Mult=IntDim                    ! Multiplicity of integral.
      N=MAX(2,Mult)                  ! "Pseudo-multiplicity" at Mult=1.
      Jump=MIN(2,Mult)               ! Switch 'Mult = 1 <--> Mult > 1'.
      N0=2*N+3                       !
      N1=N0-1                        !
      N2=N0-2                        ! Ancillari constants.
      N0dupl=2*N0                    !
      Cardin=2**N                    !
      V1=(729+(50*N-950)*N)*C1       ! V1,...,V4  and  W1,...,W5 are the
      V3=(53-20*N)*C3                ! weights for two integration rules
      W1=(3206+(100*N-2280)*N)*C0    ! B'F and BF of degree of exactness
      W3=(455-100*N)*C0              ! 5 and 7, respectively, normalized
      W5=C5/Cardin                   ! to the factor 2**N (see Ref.[1]).

      NRcall=2**N+2*N*(N+1)+1        ! Number of rule B'F calls.
      MinMax=NRcall &                !
           * INT(2*One*Length/N0-1)  ! MIN(MaxPTS) for the given Length.

      MinPTS=MinCal                  ! The minimum number of  Fun calls;
                                     ! simultaneously, MinPTS prescribes
      MaxPTS=MinMax+MinPTS           ! the maximum number of  Fun calls.

      MinFin=MaxPTS                  ! Initial value for 'MulInf'      ! 
      MaxFin=MinPTS                  ! Initial value for 'MulInf'      ! 

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!                                                                      c
!  If parameters MaxPTS or Length proved to be inadequate to the req-  c
!  uested accuracy RelErr, the routine prints an error message and is  c
!  interrupted with alternative RETURN.                                c
!                                                                      c
!  MinFin and MaxFin contain on exit the minimum and maximum numbers   c
!  of integrand evaluations performed within a given series of integ-  c
!  rations. To print these numbers it will suffice to insert operator  c
!  'CALL MULINF' in any point after that but before the next series.   c
!                                                                      c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      IF(Mult.EQ.1) NRcall=0
      Lim1=MaxPTS-2*NRcall
      Lim2=Length-N0
      RETURN

!==================================================================
      ENTRY MulTot(Fun,Res,RelErr,Aerr,MinCal)
!==================================================================

      EPS=RelErr
      Lim1=Lim1-MinPTS+MinCal
      MinPTS=MinCal
      MaxPTS=MinMax+MinPTS

!==================================================================
      ENTRY MulInt(Fun,Res)
!==================================================================

      Flag=.FALSE.
      NFcall=0
      NsubR=N0
      LsubR=N0
      Output=Zero
      AbsErr=Zero
      ALLOCATE(STORE(Length))
      IF(Mult.GT.1) GOTO 2
      CTR1=(Xupp(1)+Xlow(1))*Half
      WTH1= Xupp(1)-CTR1
      WTH2= Half
    1 VolRGN=4*WTH1*WTH2
      X(1)=CTR1
      Sum1=Fun(X)
       V=WTH1*Y2
      X(1)=CTR1-V
      Sum2=Fun(X)
      X(1)=CTR1+V
      Sum2=Sum2+Fun(X)
       V=WTH1*Y4
      X(1)=CTR1-V
      Sum4=Fun(X)
      X(1)=CTR1+V
      Sum4=Sum4+Fun(X)
      
      IF(7*Sum2.EQ.Sum4+12*Sum1) THEN
         Ndiv=2
      ELSE
         Ndiv=1
      ENDIF
      
      V=2*Sum1
      Sum2=V+Sum2
      Sum3=V+Sum4
      Sum4=2*Sum4
      
      IF(WTH1.EQ.Zero) THEN
         NFcall=NFcall+5
         IF(WTH2.GT.Zero) THEN
            Sum5=Sum1
         ELSE
            Sum5=V
         ENDIF
         GOTO 13
      ENDIF

      V=WTH1*Y5
      X(1)=CTR1-V
      Sum5=Fun(X)
         
      IF(WTH1.GT.Zero) THEN
         X(1)=CTR1+V
         NFcall=NFcall+7
         IF(WTH2.GT.Zero) THEN
            Sum5=2*(Sum5+Fun(X))
         ELSE
            Sum5=Sum5+Fun(X)
         ENDIF
      ELSE
         IF(WTH2.GT.Zero) THEN
            X(1)=CTR1+V
            Sum5=2*Sum5+Fun(X)
            NFcall=NFcall+7
         ELSE
            NFcall=NFcall+6
         ENDIF
      ENDIF

      GOTO 13

    2 CONTINUE
      DO k=1,N
         CTR(k)=(Xupp(k)+Xlow(k))*Half
         WTH(k)= Xupp(k)-CTR(k)
      ENDDO

    4 VolRGN=Cardin
    
      DO k=1,N
         VolRGN=VolRGN*WTH(k)
         X(k)=CTR(k)
      ENDDO

      Sum1=Fun(X)
      Sum2=Zero
      Sum3=Zero
      DifMax=Zero
      DO k=1,N
         V=CTR(k)
         VTH(k)=WTH(k)*Y2
         X(k)=V-VTH(k)
         S2=Fun(X)
         X(k)=V+VTH(k)
         S2=S2+Fun(X)
         Sum2=Sum2+S2
         VTH(k)=WTH(k)*Y4
         X(k)=V-VTH(k)
         S3=Fun(X)
         X(k)=V+VTH(k)
         S3=S3+Fun(X)
         Sum3=Sum3+S3
         Dif=ABS(7*S2-S3-12*Sum1)
         DifMax=MAX(Dif,DifMax)
         IF(DifMax.EQ.Dif) Ndiv=k
         X(k)=V
      ENDDO

      Sum4=Zero

      DO j0=2,N
         j=j0-1
         DO k=j0,N
            DO j1=1,2
               VTH(j)=    -VTH(j)
               X(j)=CTR(j)+VTH(j)
               DO k1=1,2
                  VTH(k)=    -VTH(k)
                  X(k)=CTR(k)+VTH(k)
                  Sum4=Sum4+Fun(X)
               ENDDO
            ENDDO
            X(k)=CTR(k)
         ENDDO
         X(j)=CTR(j)
      ENDDO

      Sum5=Zero
      DO j=1,N
         VTH(j)=    -WTH(j)*Y5
         X(j)=CTR(j)+VTH(j)
      ENDDO

   11 Sum5=Sum5+Fun(X)
      DO j=1,N
         VTH(j)=    -VTH(j)
         X(j)=CTR(j)+VTH(j)
         IF(VTH(j).GT.Zero) GOTO 11
      ENDDO

      NFcall=NFcall+NRcall
   13 ErrRGN=VolRGN*(V1*Sum1+V2*Sum2+V3*Sum3+V4*Sum4)
      ValRGN=VolRGN*(W1*Sum1+W2*Sum2+W3*Sum3+W4*Sum4+W5*Sum5)
      ErrRGN=ABS(ErrRGN-ValRGN)
      Output=    Output+ValRGN
      AbsErr=    AbsErr+ErrRGN
      
      IF(Flag) THEN
   14    IsubR=2*NsubR
         IF(IsubR.GT.LsubR) GOTO 16
         IF(IsubR.LT.LsubR) THEN
            I=IsubR+N0
            IF(STORE(IsubR).LT.STORE(I)) IsubR=I
         ENDIF
         IF(ErrRGN.GE.STORE(IsubR)) GOTO 16
         DO k=0,N1
            STORE(NsubR-k)=STORE(IsubR-k)
         ENDDO
         NsubR=IsubR
         GOTO 14
      ENDIF

   15 IsubR=(NsubR/N0dupl)*N0
      IF((IsubR.GE.N0).AND.(ErrRGN.GT.STORE(IsubR))) THEN
         DO k=0,N1
            STORE(NsubR-k)=STORE(IsubR-k)
         ENDDO
         NsubR=IsubR
         GOTO 15
      ENDIF

   16 I=NsubR-2
      DO j=1,N
         IsubR=I-(j+j)
         STORE(IsubR  )=WTH(j)
         STORE(IsubR+1)=CTR(j)
      ENDDO
      STORE(NsubR  )=ErrRGN
      STORE(NsubR-1)=ValRGN
      STORE(I)      =Ndiv

      IF(Flag) THEN
         Flag=.FALSE.
         CTR(Nd)=CTR(Nd)+2*WTH(Nd)
         LsubR=LsubR+N0
         NsubR=LsubR
         GOTO (1,4) Jump
      ENDIF
      
     IF((AbsErr.LE.(EPS*ABS(Output)+AErrr)) & 
                        .AND.(NFcall.GE.MinPTS)) THEN
         DEALLOCATE(STORE)
         MinFin=MIN(MinFin,NFcall)                                      ! *
         MaxFin=MAX(MaxFin,NFcall)                                      ! *
         Res=Output                           ! Tutto va bene
         RETURN
      ENDIF
      
      IF(NFcall.GT.Lim1) THEN
         WRITE(*   ,101) Mult,MaxPTS          ! Error message (screen)   *
         WRITE(Nlog,101) Mult,MaxPTS          ! Error message (LogFile)  *
         GOTO 18
      ENDIF
      
      IF(LsubR.GT.Lim2) THEN
         WRITE(*   ,102) Mult,Length          ! Error message (screen)   *
         WRITE(Nlog,102) Mult,Length          ! Error message (LogFile)  *
         GOTO 18
      ENDIF

      Flag=.TRUE.
      NsubR=N0
      DO j=1,N
         IsubR=N2-(j+j)
         WTH(j)=    STORE(IsubR  )
         CTR(j)=    STORE(IsubR+1)
      ENDDO

      AbsErr=AbsErr-STORE(N0)
      Output=Output-STORE(N1)
      Nd    =       STORE(N2)
      WTH(Nd)=   Half*WTH(Nd)
      CTR(Nd)=CTR(Nd)-WTH(Nd)

      GOTO (1,4) Jump
   18 CONTINUE

      DEALLOCATE(STORE) 

      IF(Output.NE.Zero) THEN                                           ! *
         WRITE(Nlog,103) AbsErr/ABS(Output),NFcall                      ! *
      ELSE                                                              ! *
         WRITE(Nlog,104) AbsErr,NFcall                                  ! *
      ENDIF                                                             ! *
      RETURN

!==================================================================
      ENTRY MulInf
!==================================================================

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!                                                                      c
!  This entry together with all starred operators (see above)  may be  c
!  removed if the output (see FORMATs) is not warranted. However that  c
!  information is helpful in the debug stage, specifically for reduc-  c
!  tion of CPU time and volume occupied in the main memory.            c
!                                                                      c
!           An example of output (format is changed a bit)             c
!                                                                      c
!  ==================================================================  c
!  |           MULINFORM (multiplicity of integral is  2)           |  c
!  |----------------------------------------------------------------|  c
!  | RelErr = 1.000D-05 | MaxCal =     9714359 | MinCal =       100 |  c
!  | Length =   2000000 | MaxFin =        6341 | MinFin =      1377 |  c
!  ==================================================================  c
!                                                                      c
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      WRITE(*,105) Mult,EPS,MaxPTS,MinPTS,Length,MaxFin,MinFin
      WRITE(*,'(A8,E12.5)')  ' Aerrr = ', Aerrr
      WRITE(Nlog,105) Mult,EPS,MaxPTS,MinPTS,Length,MaxFin,MinFin
         MinFin=MaxPTS
         MaxFin=MinPTS
      RETURN

  101 FORMAT(/ '  MULT =',I3/'  THE MAXIMUM NUMBER OF INTEGRAND CALLS,' &
          , I15,', IS TOO SMALL'/'  FOR THE REQUIRED ACCURACY:')
  102 FORMAT(/ '  MULT =',I3/'  THE STORAGE ARRAY LENGTH,'              &
          , I10,', IS TOO SMALL'/'  FOR THE REQUIRED ACCURACY:')  
  103 FORMAT('  RELATIVE ERROR =',1pd9.2                                &
                                , ', NUMBER OF INTEGRAND CALLS =',I14/)
  104 FORMAT('  RESULT = 0! ABSOLUTE ERROR =',1pd9.2                    &
                                , ', NUMBER OF INTEGRAND CALLS =',I14/)
  105 FORMAT(/16x,'MULINFORM (multiplicity of integral is',I4,')'/      &
               3x,69('-')/' | RelErr =',1PD10.3,' | MaxCal =',I15,      &
                          ' | MinCal =',    I12,' |'/                   &
                          ' | Length =',    I10,' | MaxFin =',I15,      &
                          ' | MinFin =',    I12,' |'/3x,69('-')/)

      END SUBROUTINE MulSet
!==================================================================================


!---------------------------
end module gazlib

